{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cd13daf5",
   "metadata": {},
   "source": [
    "# Parallel Thermosyphon-Driven flow\n",
    "\n",
    "\n",
    "In this example we'll demonstrate how to build parallel hydraulic loops driven by temprature differences in each vertical branch, which induce thermosyphonic flow. \n",
    "For the sake of simplicity we'll consider only the Kirchhoff equations without energy conservation.\n",
    "\n",
    "Assuming the case of three hydraulic resistors connected in parallel in vertical channels. each with different temprature (hot, mean and cold). \n",
    "This serves as a test case for inner-loop circulations that can happen in thermosyphon driven flow in reactor cores with channels at different temperatures.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\dot{m}_c&=\\dot{m}_h+\\dot{m}_m \\\\\n",
    "\\dot{m}_cK_c+\\dot{m}_hK_h&=(\\rho_c-\\rho_h)g\\Delta h \\\\\n",
    "\\dot{m}_cK_c+\\dot{m}_mK_m&=(\\rho_c-\\rho_m)g\\Delta h\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "which can be simplified into:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix}\n",
    "K_c+K_h & K_c \\\\\n",
    "K_c &K_c+K_m\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\dot{m}_h \\\\\n",
    "\\dot{m}_m\n",
    "\\end{pmatrix}=g\\Delta h\n",
    "\\begin{pmatrix}\n",
    "\\rho_c-\\rho_h \\\\\n",
    "\\rho_c-\\rho_m\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8cc7070-e9fa-4189-b3f6-11f5a58b2ed8",
   "metadata": {},
   "source": [
    "All that's left is to invert the matrix and we'll obtain the mass flow rates through each hydraulic resistor.\n",
    "This can be done by hand.\n",
    "\n",
    "The solution is:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\dot{m}_h &= \\Delta hg\\frac{(K_c+K_m)(\\rho_c-\\rho_h) - K_c(\\rho_c-\\rho_m)}{K_cK_h+K_cK_m+K_hK_m} \\\\\n",
    "\\dot{m}_m &= \\Delta hg\\frac{(K_c+K_h)(\\rho_c-\\rho_m) - K_c(\\rho_c-\\rho_h)}{K_cK_h+K_cK_m+K_hK_m} \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c778f21-00b0-483a-98db-a24e976b04bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytic_solution(dh, g, kc, kh, km, rc, rh, rm):\n",
    "    denominator = kc * kh + kc * km + kh * km\n",
    "    mh = dh * g * ((kc + km) * (rc - rh) - kc * (rc - rm)) / denominator\n",
    "    mm = dh * g * ((kc + kh) * (rc - rm) - kc * (rc - rh)) / denominator\n",
    "    return mh, mm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96f6369e",
   "metadata": {},
   "source": [
    "After we obtained the solution analytiaclly, we'll solve this problem using the `stream` library.\n",
    "\n",
    "The first step is to import the relevant objects, representing hydraulic resistors, gravity, junctions, flow graph etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d9c33b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from stream.calculations import Gravity, Resistor\n",
    "from stream.calculations.kirchhoff import Junction\n",
    "from stream.composition import flow_edge, flow_graph, flow_graph_to_agr_and_k\n",
    "from stream.substances import light_water as h2o\n",
    "from stream.units import g"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c00348c-c6f5-45f0-8b7a-86f9a1d7aa75",
   "metadata": {},
   "source": [
    "We will also need to use actual numbers for the numeric solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8276b30c-12b2-4895-a7fc-ce89f13a28fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "Tc, Th, Tm = 20.0, 50.0, 40.0\n",
    "rc, rh, rm = [h2o.density(t) for t in [Tc, Th, Tm]]\n",
    "kc, kh, km = 1, 10, 10\n",
    "dh = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "595fcb36",
   "metadata": {},
   "source": [
    "In our problem we need to create 2 junction objects, which we name \"top\" and \"bottom\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e27b69e",
   "metadata": {},
   "outputs": [],
   "source": [
    "top = Junction(\"top\")\n",
    "bottom = Junction(\"bottom\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71268fba",
   "metadata": {},
   "source": [
    "for functions and classes that require the same repeated inputs we can use `functools.partial()` function to create a partially applied function that is imbued with the repeating inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db5a9236",
   "metadata": {},
   "outputs": [],
   "source": [
    "gravity = partial(Gravity, fluid=h2o, gravity=g)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ca918d3",
   "metadata": {},
   "source": [
    "Now we can create the hydraulic circuit using `flow_graph` where each branch is defined using `flow_edge`.\n",
    "\n",
    "This is the basic way to represent hydraulic relations in `Stream`. It's important to notice that the hyraulic circuit is closed, so Kirchhoff calculatiosn can do closed loop equations. Hence, for the example below, if you defined `flow_edge` that started from a `junction` named `top` you have to define a `flow_edge` returning to `top`. Another important point here is the sign of `Gravity`'s object disposition $\\Delta h$. In `Stream` it should be positive for downard flow (in the direction of gravity) and negative for upward flow. The flow rates are positive when the flow is in the same direction as the appropriate flow edge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b181f0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fg = flow_graph(\n",
    "    flow_edge((top, bottom), Resistor(kh, \"Rh\"), gravity(disposition=dh, name=\"gh\")),\n",
    "    flow_edge((top, bottom), Resistor(km, \"Rm\"), gravity(disposition=dh, name=\"gm\")),\n",
    "    flow_edge((bottom, top), Resistor(kc, \"Rc\"), gravity(disposition=-dh, name=\"gc\")),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "452be2e9",
   "metadata": {},
   "source": [
    "The next step is to build an aggregator of the equations out of the `flow_graph` we constructed, which can be done using `flow_graph_to_agr_k`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "838e917f",
   "metadata": {},
   "outputs": [],
   "source": [
    "agr, _ = flow_graph_to_agr_and_k(fg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc1cd9c1",
   "metadata": {},
   "source": [
    "All that's left before solving the problem is to define the temperatures at each branch. \n",
    "Because we only look at momentum conservation and not really intrested in energy conservation in this toy problem, it's enough to define the temperatures of each gravity object using the concept of `ext_funcs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5d8edd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "agr.funcs = {\n",
    "    agr[\"gh\"]: {\"Tin\": Th, \"Tin_minus\": Th},\n",
    "    agr[\"gm\"]: {\"Tin\": Tm, \"Tin_minus\": Tm},\n",
    "    agr[\"gc\"]: {\"Tin\": Tc, \"Tin_minus\": Tc},\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a625e42c",
   "metadata": {},
   "source": [
    "A very nice tool to inspect the number of equations and variables is the `report` function. This will tell us if something is not properlly defined or missing before trying to solve the agregator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f01f35f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from stream.analysis.report import report\n",
    "\n",
    "report(agr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c42fd257",
   "metadata": {},
   "source": [
    "Finally we'll solve the problem using the `agr.solve_steady()` method with an initial guess of a vector of ones. \n",
    "We also use `agr.save()` to transform the solution from a vector to a dictionary with keys that indicate the context of each value in the result vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21133646",
   "metadata": {},
   "outputs": [],
   "source": [
    "stream_sol = agr.save(agr.solve_steady(np.ones(len(agr)), options={\"xtol\": 1e-13}))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3f1d48e-3215-428c-8f12-479fc1dab0d4",
   "metadata": {},
   "source": [
    "The analytic function should also be computed for our specific values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1d44072-03b0-4734-94c2-3e5c7068cabd",
   "metadata": {},
   "outputs": [],
   "source": [
    "analytic_sol = analytic_solution(dh, g, kc, kh, km, rc, rh, rm)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1a5c379-29e2-4523-9c5c-c8511ace2ab7",
   "metadata": {},
   "source": [
    "Then we can show what the solution is and `stream`'s error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb4ad738",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdot_h_stream = stream_sol[\"Kirchhoff\"][\"(top -> bottom, 0)\"]\n",
    "mdot_h_analytical = -analytic_sol[0]\n",
    "print(f\"mdot_h={mdot_h_analytical:.3e}\\nerror:{100 * (mdot_h_analytical - mdot_h_stream) / mdot_h_analytical:.5e} %\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a701fffa",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdot_m_stream = stream_sol[\"Kirchhoff\"][\"(top -> bottom, 1)\"]\n",
    "mdot_m_analytical = -analytic_sol[1]\n",
    "print(f\"mdot_m={mdot_m_analytical:.3e}\\nerror:{100 * (mdot_m_analytical - mdot_m_stream) / mdot_m_analytical:.5e} %\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29f5ddc9",
   "metadata": {},
   "source": [
    " Thus the numeric result obtained is really close (numerically identical) to the analytical solution"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}